function ng = pij (i,j)
  global N
  if(i < 4*N + 1) ng = j - N + N*( i - 1);
  else ng = j + 2*N*(i - 1) - 4*N^2;
  end
end

